!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine regrid(current)
!==================================================================================================================================|

   use parameters
   use block_define
   use interface_block_generate
   use interface_hilbert_curve
   implicit none

   type(block), pointer, intent(inout) :: current

   integer(4) :: i_block, j_block, k_block, i_neigh, j_neigh, k_neigh, ip1, jp1, kp1, im1, jm1, km1
   integer(4) :: flag, flag_destroy, n_blocks_tmp, n_gene_tmp
   integer(4), dimension(:, :, :), allocatable :: gene_tmp

   real(8) :: lv_exp, num_tmp, thral1
   real(8), dimension(n_thr) :: dva, dva2, d2va2
   real(8), dimension(nx_block, ny_block, nz_block) :: vx, vy, vz, pr, te

   type(block), pointer :: child, neigh, next, framework_next, pointer_tmp

!----------------------------------------------------------------------------------------------------------------------------------|
!  error estimation (flag = 1 then block generation, flag = -1 then block destruction and flag = 0 for nothing)
!----------------------------------------------------------------------------------------------------------------------------------|
   flag = 0

   if (current%lv .lt. n_levels) then

      flag = -1

      lv_exp = 2.0d0 ** (thr_ex * (current%lv - 1))
      num_tmp = sqrt(3.0d0) / 2.0d0
      thral1 = 1.0d0 - thr_al
      call variables_transform (vx, vy, vz, pr, te,                                                                                &
         current%ro, current%mx, current%my, current%mz, current%bx, current%by, current%bz, current%en,                           &
         nx_block, ny_block, nz_block, 1)
      loop3: do k = 2, nz_block - 1
         kp1 = k + 1
         km1 = k - 1
         do j = 2, ny_block - 1
            jp1 = j + 1
            jm1 = j - 1
            do i = 2, nx_block - 1
               ip1 = i + 1
               im1 = i - 1

               dva2(1) = num_tmp * sqrt((current%ro(ip1, j, k) - current%ro(im1, j, k)) ** 2 +                                     &
                                        (current%ro(i, jp1, k) - current%ro(i, jm1, k)) ** 2 +                                     &
                                        (current%ro(i, j, kp1) - current%ro(i, j, km1)) ** 2)
               dva2(2) = num_tmp * sqrt((vx(ip1, j, k) - vx(im1, j, k)) ** 2 +                                                     &
                                        (vx(i, jp1, k) - vx(i, jm1, k)) ** 2 +                                                     &
                                        (vx(i, j, kp1) - vx(i, j, km1)) ** 2)
               dva2(3) = num_tmp * sqrt((vy(ip1, j, k) - vy(im1, j, k)) ** 2 +                                                     &
                                        (vy(i, jp1, k) - vy(i, jm1, k)) ** 2 +                                                     &
                                        (vy(i, j, kp1) - vy(i, j, km1)) ** 2)
               dva2(4) = num_tmp * sqrt((vz(ip1, j, k) - vz(im1, j, k)) ** 2 +                                                     &
                                        (vz(i, jp1, k) - vz(i, jm1, k)) ** 2 +                                                     &
                                        (vz(i, j, kp1) - vz(i, j, km1)) ** 2)
               dva2(5) = num_tmp * sqrt((current%bx(ip1, j, k) - current%bx(im1, j, k)) ** 2 +                                     &
                                        (current%bx(i, jp1, k) - current%bx(i, jm1, k)) ** 2 +                                     &
                                        (current%bx(i, j, kp1) - current%bx(i, j, km1)) ** 2)
               dva2(6) = num_tmp * sqrt((current%by(ip1, j, k) - current%by(im1, j, k)) ** 2 +                                     &
                                        (current%by(i, jp1, k) - current%by(i, jm1, k)) ** 2 +                                     &
                                        (current%by(i, j, kp1) - current%by(i, j, km1)) ** 2)
               dva2(7) = num_tmp * sqrt((current%bz(ip1, j, k) - current%bz(im1, j, k)) ** 2 +                                     &
                                        (current%bz(i, jp1, k) - current%bz(i, jm1, k)) ** 2 +                                     &
                                        (current%bz(i, j, kp1) - current%bz(i, j, km1)) ** 2)
               dva2(8) = num_tmp * sqrt((pr(ip1, j, k) - pr(im1, j, k)) ** 2 +                                                     &
                                        (pr(i, jp1, k) - pr(i, jm1, k)) ** 2 +                                                     &
                                        (pr(i, j, kp1) - pr(i, j, km1)) ** 2)
               dva2(9) = num_tmp * sqrt((te(ip1, j, k) - te(im1, j, k)) ** 2 +                                                     &
                                        (te(i, jp1, k) - te(i, jm1, k)) ** 2 +                                                     &
                                        (te(i, j, kp1) - te(i, j, km1)) ** 2)

               d2va2(1) = sqrt((current%ro(ip1, j, k) - 2.0d0 * current%ro(i, j, k) + current%ro(im1, j, k)) ** 2 +                &
                               (current%ro(i, jp1, k) - 2.0d0 * current%ro(i, j, k) + current%ro(i, jm1, k)) ** 2 +                &
                               (current%ro(i, j, kp1) - 2.0d0 * current%ro(i, j, k) + current%ro(i, j, km1)) ** 2)
               d2va2(2) = sqrt((vx(ip1, j, k) - 2.0d0 * vx(i, j, k) + vx(im1, j, k)) ** 2 +                                        &
                               (vx(i, jp1, k) - 2.0d0 * vx(i, j, k) + vx(i, jm1, k)) ** 2 +                                        &
                               (vx(i, j, kp1) - 2.0d0 * vx(i, j, k) + vx(i, j, km1)) ** 2)
               d2va2(3) = sqrt((vy(ip1, j, k) - 2.0d0 * vy(i, j, k) + vy(im1, j, k)) ** 2 +                                        &
                               (vy(i, jp1, k) - 2.0d0 * vy(i, j, k) + vy(i, jm1, k)) ** 2 +                                        &
                               (vy(i, j, kp1) - 2.0d0 * vy(i, j, k) + vy(i, j, km1)) ** 2)
               d2va2(4) = sqrt((vz(ip1, j, k) - 2.0d0 * vz(i, j, k) + vz(im1, j, k)) ** 2 +                                        &
                               (vz(i, jp1, k) - 2.0d0 * vz(i, j, k) + vz(i, jm1, k)) ** 2 +                                        &
                               (vz(i, j, kp1) - 2.0d0 * vz(i, j, k) + vz(i, j, km1)) ** 2)
               d2va2(5) = sqrt((current%bx(ip1, j, k) - 2.0d0 * current%bx(i, j, k) + current%bx(im1, j, k)) ** 2 +                &
                               (current%bx(i, jp1, k) - 2.0d0 * current%bx(i, j, k) + current%bx(i, jm1, k)) ** 2 +                &
                               (current%bx(i, j, kp1) - 2.0d0 * current%bx(i, j, k) + current%bx(i, j, km1)) ** 2)
               d2va2(6) = sqrt((current%by(ip1, j, k) - 2.0d0 * current%by(i, j, k) + current%by(im1, j, k)) ** 2 +                &
                               (current%by(i, jp1, k) - 2.0d0 * current%by(i, j, k) + current%by(i, jm1, k)) ** 2 +                &
                               (current%by(i, j, kp1) - 2.0d0 * current%by(i, j, k) + current%by(i, j, km1)) ** 2)
               d2va2(7) = sqrt((current%bz(ip1, j, k) - 2.0d0 * current%bz(i, j, k) + current%bz(im1, j, k)) ** 2 +                &
                               (current%bz(i, jp1, k) - 2.0d0 * current%bz(i, j, k) + current%bz(i, jm1, k)) ** 2 +                &
                               (current%bz(i, j, kp1) - 2.0d0 * current%bz(i, j, k) + current%bz(i, j, km1)) ** 2)
               d2va2(8) = sqrt((pr(ip1, j, k) - 2.0d0 * pr(i, j, k) + pr(im1, j, k)) ** 2 +                                        &
                               (pr(i, jp1, k) - 2.0d0 * pr(i, j, k) + pr(i, jm1, k)) ** 2 +                                        &
                               (pr(i, j, kp1) - 2.0d0 * pr(i, j, k) + pr(i, j, km1)) ** 2)
               d2va2(9) = sqrt((te(ip1, j, k) - 2.0d0 * te(i, j, k) + te(im1, j, k)) ** 2 +                                        &
                               (te(i, jp1, k) - 2.0d0 * te(i, j, k) + te(i, jm1, k)) ** 2 +                                        &
                               (te(i, j, kp1) - 2.0d0 * te(i, j, k) + te(i, j, km1)) ** 2)

               dva(1) = (thr_al * dva2(1) / (abs(current%ro(i, j, k)) + eps) +                                                     &
                         thral1 * d2va2(1) / (dva2(1) + thr_fl * (abs(current%ro(i, j, k)) + eps))) * lv_exp
               dva(2) = (thr_al * dva2(2) / (abs(vx(i, j, k)) + 0.1d0) +                                                           &
                         thral1 * d2va2(2) / (dva2(2) + thr_fl * (abs(vx(i, j, k)) + 0.1d0))) * lv_exp
               dva(3) = (thr_al * dva2(3) / (abs(vy(i, j, k)) + 0.1d0) +                                                           &
                         thral1 * d2va2(3) / (dva2(3) + thr_fl * (abs(vy(i, j, k)) + 0.1d0))) * lv_exp
               dva(4) = (thr_al * dva2(4) / (abs(vz(i, j, k)) + 0.1d0) +                                                           &
                         thral1 * d2va2(4) / (dva2(4) + thr_fl * (abs(vz(i, j, k)) + 0.1d0))) * lv_exp
               dva(5) = (thr_al * dva2(5) / (abs(current%bx(i, j, k)) + 0.1d0) +                                                   &
                         thral1 * d2va2(5) / (dva2(5) + thr_fl * (abs(current%bx(i, j, k)) + 0.1d0))) * lv_exp
               dva(6) = (thr_al * dva2(6) / (abs(current%by(i, j, k)) + 0.1d0) +                                                   &
                         thral1 * d2va2(6) / (dva2(6) + thr_fl * (abs(current%by(i, j, k)) + 0.1d0))) * lv_exp
               dva(7) = (thr_al * dva2(7) / (abs(current%bz(i, j, k)) + 0.1d0) +                                                   &
                         thral1 * d2va2(7) / (dva2(7) + thr_fl * (abs(current%bz(i, j, k)) + 0.1d0))) * lv_exp
               dva(8) = (thr_al * dva2(8) / (abs(pr(i, j, k)) + eps) +                                                             &
                         thral1 * d2va2(8) / (dva2(8) + thr_fl * (abs(pr(i, j, k)) + eps))) * lv_exp
               dva(9) = (thr_al * dva2(9) / (abs(te(i, j, k)) + eps) +                                                             &
                         thral1 * d2va2(9) / (dva2(9) + thr_fl * (abs(te(i, j, k)) + eps))) * lv_exp

               if (any(dva .gt. thr_va0)) then
                  flag = 1
                  exit loop3
               else if (any(dva .gt. thr_va0 * thr_de)) then
                  flag = 0
               endif
            enddo
         enddo 
      enddo loop3
   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  grid generation or destruction
!----------------------------------------------------------------------------------------------------------------------------------|
   n_blocks_tmp = n_blocks(cpu_id_p1)
   if (flag .eq. 1) then

      if (current%nc .eq. 0 .and. current%lv .lt. n_levels) then
         i_gene(cpu_id_p1) = i_gene(cpu_id_p1) + 1
!----------------------------------------------------------------------------------------------------------------------------------|
!        for the case i_gene(cpu_id_p1) is greater than the n_gene
!----------------------------------------------------------------------------------------------------------------------------------|
         if (i_gene(cpu_id_p1) .gt. n_gene) then
            allocate(gene_tmp(n_dims, n_levels, n_gene))
            gene_tmp = gene
            n_gene_tmp = n_gene
            n_gene = n_gene + ini_block_x * ini_block_y * ini_block_z
            deallocate(gene)
            allocate(gene(n_dims, n_levels, n_gene))
            gene(:, :, 1 : n_gene_tmp) = gene_tmp
            deallocate(gene_tmp)
         endif
         gene(:, :, i_gene(cpu_id_p1)) = current%id
         call block_generate(current, current%cpu_id)
         next => current%next
         framework_next => current%framework_next
         pointer_tmp => current
         call hilbert_curve(current, pointer_tmp, current%parent_cell, current%hilbert_pos, 1, (/1, 1, 1/), current%cpu_id)
         next%last => pointer_tmp
         pointer_tmp%next => next
         framework_next%framework_last => pointer_tmp
         pointer_tmp%framework_next => framework_next
         max_level = max(max_level, current%lv + 1)
      endif

   elseif (flag .eq. -1) then

      flag_destroy = 1
      do k_block = 1, current%nc_z
      do j_block = 1, current%nc_y
      do i_block = 1, current%nc_x
         child => current%child(i_block, j_block, k_block)%pblock
         do k_neigh = -1, 1
         do j_neigh = -1, 1
         do i_neigh = -1, 1
            neigh => child%neigh(i_neigh, j_neigh, k_neigh)%pblock
            if (associated(neigh)) then
               if (all(neigh%parent%id .eq. child%parent%id)) cycle
               if (neigh%lv .eq. child%lv .and. neigh%nc .ne. 0) flag_destroy = 0
            endif
         enddo
         enddo
         enddo
      enddo
      enddo
      enddo

      if (flag_destroy .eq. 1 .and. current%nc .ne. 0) then
         i_gene(cpu_id_p1) = i_gene(cpu_id_p1) + 1
!----------------------------------------------------------------------------------------------------------------------------------|
!        for the case i_gene(cpu_id_p1) is greater than the n_gene
!----------------------------------------------------------------------------------------------------------------------------------|
         if (i_gene(cpu_id_p1) .gt. n_gene) then
            allocate(gene_tmp(n_dims, n_levels, n_gene))
            gene_tmp = gene
            n_gene_tmp = n_gene
            n_gene = n_gene + ini_block_x * ini_block_y * ini_block_z
            deallocate(gene)
            allocate(gene(n_dims, n_levels, n_gene))
            gene(:, :, 1 : n_gene_tmp) = gene_tmp
            deallocate(gene_tmp)
         endif
         gene(:, :, i_gene(cpu_id_p1)) = - current%id
      endif

   endif

!----------------------------------------------------------------------------------------------------------------------------------|
!  reset n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
   n_blocks(cpu_id_p1) = n_blocks_tmp

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine regrid
